Exact solution of the ID Hubbard model in the atomic limit 
with inter-site magnetic coupling 
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In this paper we present for the first time the exact solution in the narrow-band limit of the ID 
extended Hubbard model with nearest-neighbour spin-spin interactions described by an exchange 
constant J. An external magnetic field h is also taken into account. This result has been obtained in 
the framework of the Green's functions formalism, using the Composite Operator Method. By means 
of this theoretical background, we have studied some relevant features such as double occupancy, 
magnetization, spin-spin and charge-charge correlation functions and derived a phase diagram for 
both ferro (J > 0) and anti-ferro (J < 0) coupling in the limit of zero temperature. We also 
report a study on density of states, specific heat, charge and spin susceptibilities. In the limit of 
zero temperature, we show that the model exhibits a very rich phase diagram characterized by 
different magnetic orders and by the coexistence of charge and spin orderings at commensurate 
filling. Moreover, our analysis at finite temperature of density of states and response functions 
shows the presence of low-temperature charge and spin excitations near the phase boundaries. 
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PACS numbers: 71.27.+a,75.10.Pq,71.10.Hf 

I. INTRODUCTION 

For many years the Mott-Hubbard theory has been 
considered as the archetype for a class of Hamiltonian 
models aiming the description of correlations between 
electrons in the d-bands of transition metal compounds. 
However, although the Hubbard model^ and its reduc- 
tion, the t-J model 2 , have been successful in repro- 
ducing a plethora of different anomalous phenomena 
among which magnetic orders^—, Metal Insulator Tran- 
sition (MIT) 7 ., Spin Density Waves (SDW) 8 and High- 
Temperature Superconductivity (HTSC)^r— , their appli- 
cation to a large class of compounds is still controversial. 
The Mott-Hubbard theory in its simplest formulation, 
leads to a band gap of the order U (~ 7-10eV in oxides) 
which is difficult to justify for charge-transfer insulators 
such as Co, Ni and Cu. Furthermore, it is also difficult to 
understand the mctallicity of many sulfides among which 
NiS, CuS and CoS which would require a strong reduc- 
tion of the on-site interaction U to l-2eV— . 

These and many other issues led to the conviction that 
single-band Hubbard and t-J models are not sufficient 
to catch all the relevant features due to electronic cor- 
relations. Therefore, with the aim of describing strong 
electron-electron correlations as well as interactions of 
electrons with other degrees of freedoms (such as lattice 
vibrations, light), several extensions of the bare Hub- 
bard Hamiltonian have been proposed resulting in the 
introduction of the so-called Extended Hubbard Models 
(EHMs). Among these, in the recent years an increas- 
ing interest arose in the study of the t-U-J model as the 
minimal model capable to reproduce the exchange cor- 
relations, widely believed to be the basis of the pairing 
mechanism in cuprates^r— . Unlikely the Hubbard and 
t-J models, in the t-U-J Hamiltonian the exchange cou- 



pling J is not related to the Hubbard U as J » At 2 /U , 
allowing finite exchange correlations even in the pres- 
ence of strong on-site couplings. Furthermore, contrarily 
to the t-J model, an independent treatment of J does 
not necessarily require the U — > oo limit in which charge 
fluctuations are heavily suppressed. 

The t-U-J model has been analyzed in ID and 2D 
cases. In the 2D case the t-U-J Hamiltonian has been 
used as the minimal model capable to describe the 
charge-transfer nature of cuprates. In particular, sev- 
eral studies have been done to understand the influence 
of the on-site Coulomb repulsion and spin-spin exchange 
on superconductivit y 18 ' 19 . 

Motivated by the discovery of close proximity of mag- 
netic and superconducting ordering in (TMTSF^X fam- 
ily of quasi one-dimensional Bechgaard salts2S, analytical 
and numerical studies on the t-U-J model have also been 
performed in the one-dimensional case. At half-filling, 
the ground-state phase diagram of the ID t-U-J model 
has been intensively studied for both ferro and anti-ferro 
magnetic couplings in the weak-coupling limit in which 
U, J <C t. By means of the bosonization procedure, it has 
been pointed out that, in the presence of ferromagnetic 
interactions, the system is dominated by superconducting 
and spin-density-wave instabilities even in the presence 
of moderate values of the on-site Hubbard interaction^. 
On the contrary, in the presence of an anti-ferromagnetic 
exchange, bosonization and transfer-matrix renormaliza- 
tion group methods showed that at half-filling the ground 
state of the system is a Mott insulator characterized by 
spontaneous dimerization for U <C J. A transition to 
a gapless spin liquid phase occurs at U c » J/2 21 . The 
ID t-U-J model at n = 1 has also been studied with the 
inclusion of an easy-plane anisotropy in the exchange in- 
teraction in order to investigate the coexistence of triplet 
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superconductivity and ferromagnetism in a class of quasi- 
one-dimensional material a 22 ' 23 . It has been shown that, 
in the large bandwidth limit, magnetic correlations are 
enhanced by the presence of a repulsive Coulomb poten- 
tial U and a transverse spin-exchange interaction between 
electrons on nearest-neighbor (NN) sites. Therefore the 
coexistence of antiferromagnetism and triplet supercon- 
ductivity is no longer observed except for small values 
of the Coulomb interaction 2 - 2 -. Recently, it has been also 
shown that in the weak-coupling limit and for n = 1, in 
the presence of an isotropic anti-ferromagnetic exchange, 
CDW and bond SDW phases are suppressed and the 
ground state exhibits an insulating behavior character- 
ized by SDW and bond CDW phases^. 

Because of the complexity of the model, in spite of 
numerous attempts, there are no exact solutions for the 
t-U-J Hamiltonian and only few results are known out- 
side of the weak coupling regime or away from half- filling. 
Within this context we present in this work the exact so- 
lution of the ID t-U-J model in the atomic limit and in 
the presence of an external magnetic field h. It is impor- 
tant to stress that the t-U-J Hamiltonian, even in the 
atomic limit, does not represent an abstract toy model 
since a number of quasi-one dimensional insulating com- 
pounds, like CsMnBr 3 , Sr 3 ZnIrO, CuGe0 3 , Ca 3 CoRh0 6 
and SrsCuPtadri^Oe show magnetic features that can 
be described by the introduction of an inter-site mag- 
netic coupling—^. Moreover, because of charge-charge 
and spin-spin interaction terms, the U-J Hamiltonian in 
the presence of an external magnetic field represents a key 
model for the study of phase transitions between different 
magnetic orders or a good starting point for perturbativc 
approaches in terms of the hopping parameter t. Be- 
sides all the motivations given above, we would like also 
to stress that the model is interesting from the point of 
view of statistical mechanics, showing several phase tran- 
sitions, two tricritical points, anomalous behaviors for 
all the system response functions near the phase bound- 
aries. Furthermore, as we shall see in Sec HID I for the one- 
dimensional case, the U-J-h model can also be mapped 
to a spin ladder Hamiltonian with effective inter-chain 
and intra-chain spin-spin interactions, allowing for the 
exact solution of a quasi two-dimensional spin system. 

Hereafter we report a summary of the work. In the sec- 
ond Section we introduce the Green's functions formalism 
and present the exact solution of the model obtained us- 
ing the Composite Operator Metho d 27 ' 28 which, in the 
last fifteen years, has been successfully applied to sev- 
eral models and materials among which: Hubbard^ 9 - - — , 
p-cgl, i-jM, t-t'-U^, extended Hubbard (t-U-V)^., 
Kondc 37 , Anderson^, two-orbital Hubbar d 39 ' 40 , exact 
solvable models 4 ^, J 1 -J 2 ^., Hubbard-Kondo 4 ^, and 
Cuprates^ 7 -. We show that, by means of algebraic rela- 
tions, the hierarchy of the equations of motion closes and 
the Green's functions can be expressed in terms of a finite 
number of parameters to be determined self-consistently. 
A collection of results, including single particle correla- 
tion functions, charge and spin susceptibilities, density of 



states and some thermodynamic quantities, obtained in 
the limit of zero temperature and for finite T , is reported 
in the third and in the fourth Sections, respectively. 

Our study in the limit of zero temperature shows that 
the model exhibits a very rich phase diagram in which 
different magnetic orders are generated by the competi- 
tion between J, h and U in the whole range of the filling 
n. In particular, in the presence of an anti-ferromagnetic 
inter-site coupling, we observe two different types of fer- 
romagnetic order for high values of the magnetic field. 
Coexistence of charge and spin ordering has also been 
observed for commensurate fillings. Finally, our study on 
density of states at T — > and thermodynamic quantities 
at finite temperature shows clear signatures of charge and 
spin excitations in all the studied thermodynamic func- 
tions: specific heat, spin- and charge-susceptibilities. 



II. THE MODEL - COMPOSITE FIELDS AND 
GREEN'S FUNCTION 



A simple generalization of the Hubbard model can be 
obtained by including magnetic inter-site interactions. In 
this case, the Hubbard Hamiltonian can be written as 
follows: 

H = -/i%)c t (i)c(j) + ?7^n t (iK(i) + 



(1) 



where c(i) and c'(i) are annihilation and creation op- 
erators of electrons in the spinorial notation: c^(i) = 

( c|(i) c|(i) J, satisfying canonical anti-commutation 
relations. The spinorial notation will be used for all 
fermionic operators. The Heisenberg picture is used 
[i = (i, t)], i is a vector of the lattice; tij denotes the 
transfer integral and describes hopping between different 
sites; /i is the chemical potential. n a {i) = cJ.(i)c CT (i) is 
the number density operator of electrons at the site i with 
spin a. The intensity of the local Coulomb interaction is 
parameterized by U; n^(i) — n^{i) — nj.(z) is the third 
component of the spin density operator; Jj j is the ex- 
change inter-site interaction; h represents the strength of 
the external magnetic field. In this work we restrict our 
analysis to the narrow-band limit and consider only first 
neighbor interactions by taking Jj j = —2dJai_j, where d 
is the dimensionality of the system and aij is the projec- 
tion operator on the NN sites. For a d-dimensional cubic 
Bravais lattice of lattice constant a, the Fourier trans- 
form of (Xij is F.T. [ctij] — 2 Sn=l cos (k n a). Then, the 
Hamiltonian ([TJ can be written under the form: 



H = 



— nn(i) + U D(i) — hn 3 (i) — dJn 3 (i)n^ (i) 



(2) 



where n(i) — c' (i)c(i) is the total density operator 
and D(i) — n\{i)n\.{i)=hn{i)[n{i) — 1] the double oc- 
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cupancy operator. Hereafter, for a generic operator 
we use the following notation: = J^. ay$(j,t). 

We note that © is invariant under the transformation: 
(ft — >• —ft , n-f — > n^). Also, under the particle-hole trans- 
formation, the chemical potential scales as /i(2 — n) = 
U - f i(n). 

To solve the Hamiltonian ([2]), we shall use the formal- 
ism of Green's functions and equations of motion^. As 
a first step, we show that there exists a closed set of 
eigenoperators and eigenvalues of H. To this end, we 
introduce the composite field operators: 



^« = ewKwr 1 , < } « = m^d)}*- 1 , (3) 

where = [1 — n(i)]c(i) and rj(i) = n(i)c(i) are the 
Hubbard operators responsible for the transitions |0)j •O- 
\<j)i at the site i, and \a)i <H- | t4)-<i at site i, respectively. 
These fields satisfy the equations of motion: 



where cr 3 is the third Pauli matrix. By taking higher- 
order time derivatives we generate a hierarchy of com- 
posite operators. However, on the basis of the algebraic 
relations: 

n (p)(i) = n (i) + a p D(i) [n 3 {i)] 2 P = n(i) - 2D(i) 
a p = 2P-2 ' M*)] 2p+1 =»s(») 

DP{i) = D(i) n(i)n 3 (i) = n 3 (i) 

n P(i)D(i) = 2PD(i) D(i)n 3 (i) = 

for p ^ 1, the following recursion formula can be estab- 
lished for the field [nf (i)] p : 



fo) 



2p-l _ 



2^m=l ^ 



Mr 



[«§■(*) 



f(i)] a 

12m 



(6) 



where i„ are rational numbers, satisfying the sum rule 



Z^m=l Si™ 



(P) 



1, and, for 1 ^ p ^ m, A. 



(/--) 



p > m the expressions of the coefficients i„ depend on 
the coordination number z = 2d. In the Appendix [X] we 
report some values. 

We now define the composite operator: 



5 pm . For 



where and ^"(i) are multiplet operators of rank 



( ^f{i) \ 



(7) 



2z + l: 



( &•(*) \ 



77 CT (i)[nf(i)] 
1,(i)K(i)] 2 

W«K(*)] 22+1 / 



(8) 



On the basis of the equations of motion (Q} and by means 
of the recursion rules (JSJ) , it is easy to see that these fields 
are eigenoperators of the Hamiltonian ©: 



(9) 



where and are the energy matrices which can 
be calculated by means of the equations of motion (|4]) 
and the recursion rule ([5]), whose eigenvalues, Ea and 
E < £ l,n \ with n G {l,4d+ 1}, determine the quasiparticle 

excitation spectrum. Explicit expressions of and E^ 
(s = £, r/) for z = 2 are given in Appendix [Bl 

The knowledge of a complete set of eigenoperators and 
eigenvalues of the Hamiltonian allows for an exact ex- 
pression of the retarded Green's function (GF): 



G«(t-f) 



R 



V>«(M)^ s)t (M') 



(10) 



and consequently, by using the spectral theorem, of the 
correlation function (CF): 



(ii) 



In the above equations s = (£,?7), (■ • ■ ) denotes the 
quantum-statistical average over the grand canonical en- 
semble while R represents the retarded operator with 
(R[iP(x)^(y)}) = e(t x -t y ){{^(x),iPHy)}). By means 
of the field equations (J9j) , this Green's function satisfies 
the equation: 



:(') 



(s) 

where I a is the normalization matrix: 

/W = ({^)(i,t),4 s)t (i,t)} 
The solution of Eq. ifT!?)) gives: 



2z + l 



E 



(s,n) 



1 UJ — Ea +10 



(12) 



(13) 



(14) 



C«(w) = tt pi s ' n) Ti s ^S (u - E^) , (15) 



2z+l 



n=l 
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where T a s ' n) = 1 + tanh (/3E a s ' n) /2) , (3 = l/k B T and 
the spectral density matrices p a s ' n ^ are given by: 

2z+l 



(16) 



(s) 

Oct is the matrix whose columns are the eigenvectors 



of the matrix e 



Calculations of the matrices f2 



are reported in Appendix [B] It is worth noting that 

(£) (v) (s n) 

Oct = Oct • It can be shown that the pi- ' satisfy the 



sum rule £„li = 7i 

„0,n) 



The relevant elements of 



^22+1 ^(s,n) _ 

a,ab' 

Pa ' ' are given in Appendix [El for z = 2. It is important 
to stress that the expressions (|T4")) and (|15p are exact and, 
in principle, give an exact solution of the Hamiltonian ([2]) 
for any dimension of the system. However, in order to 
obtain quantitative information of the physical proper- 
ties of the system we need to have a complete knowledge 
of the Green's and correlation functions. As (fT4)) and 
(1151) show, both of them depend on two quantities: the 
energy matrix and the normalization matrix. The en- 
ergy matrices are exactly known (see Appendix [Bj . Let 
us then concentrate the attention on the normalization 
matrix. At first, we note that by using the recursion rules 
the matrix elements I^-ih can be expressed in terms 
of the elements of the first raw (p = 1, . . . , 2z + 1): 



j(£) 

r(l) 



A 



6 W - A 

(p) 



(/>) 



(p = l,...,2* + l) 



(17) 



nf(i)] 



p-i 



and A, 



where = 

n. CT (i) [n.3 (i)] p Thus, the correlation functions 

depend on the external parameters n = (n(i)), T, U, J, 
ft and on the internal parameters: /1, k^ p \ Act . It is easy 
to show that the CFs obey the following self-consistent 
equations: 

Ci% + = - A« (p = 1, . . . , 2z + 1) , (18) 



where C CT S ' ) = (ijji s \i)ipi s ^ (i)^ . The number of these 
equations is not sufficient to determine all the inter- 
nal parameters, and one needs other equations. This 
problem will be considered in the next Section, where a 
self-consistent scheme, capable of computing the internal 
parameters, will be formulated for the one-dimensional 
case. 



A. Self-consistent equations 

The previous analysis shows that the complete solution 
of the model requires the knowledge of the parameters fi, 
n( p \ and Act • These quantities may be computed by 
using algebra constraints and symmetry requirements. 



Let us consider a one-dimensional system and fix one 
site, say i, at some arbitrary point of the chain. We split 
the Hamiltonian ^ as the sum of two terms^: 



H = + H { p 
Hf ] = -2Jn 3 {i)n%{i) 



(19) 



and introduce the -representation: the statistical av- 
erage of any operator can be expressed as: 



(20) 



where (• • • )o , stands for the thermal average with re- 
spect to the reduced Hamiltonian H^: (...)„. = 



Tr{. 



. e 



}/Tr{ ( 



\. Because of transla- 



tional invariance, hereafter the dependence on the in- 
dex i will be omitted. As it is shown in Appendix [B] 
the parameters and A^ can be written as functions 
of three parameters X\ = (n a (i)) Q , X 2 = (n 3 (i)} and 
A3 = (D a (i)) Q , in terms of which one finds a solution of 
the model. In the ifo-representation, by exploiting the 
translational invariance along the chain: 



(«(*)) = (n<*(i)) 
(n 3 (i)) = <«?(*)> 
(D(i)) = (D a (i)) 



(21) 



one obtains three equations, reported in Appendix IB"! 
(|B24)) - (|B26)) . allowing to determine X\, X 2 and X 3 in 
terms of p. The chemical potential is fixed by means of 
the equation (n(i)) — n that gives: 

F 4 = (1- »)(e-^°) + G x -1 = 0, (22) 
where n is the particle number per site, and 

{e~ pH t\ = l + 2aG 2 X 2 [l + b(X 1 -2X 3 )} + 

+ (Gi - 2G 3 ) [a 2 X 2 + 2b (X x - 2X 3 ) + 
+b 2 {X 1 -2X 3 ) 2 ] . (23) 

The final expressions for G\, G 2 and G 3 are reported in 
Appendix [B] (IB161) . Thus, (|B%)) - (|B26)) and ||22J) con- 
stitute a system of coupled equations ascertaining the 
four parameters p, X\, X 2 and X 3 in terms of the ex- 
ternal parameters of the model n, U, J, h, T. Once 
these quantities are known, all local properties of the 
model can be computed. In particular, the magnetiza- 
tion m = \ (n 3 (i)) and the double occupancy D are given 
by: 



1 



2{e-f> H ? ) ) 



{X 2 + {G 1 -2G 3 )X 2 [b 



D 



+ (a 2 + b + & 2 )(Ai - 2X3)] + aG 2 [X$(l + b) + 
+ (Xt - 2A 3 ) + b{X x - 2X 3 ) 2 ]} , (24) 
G 3 



(25) 
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It is worth noting that from Eqs. (l2"TT) and (|B24I) we can 
derive the following exact relation between the particle 
number and the double occupancy 



1 = g D + n , 



(26) 



where g = (1 — G\)/G 3 . The local charge, spin and 
double occupancy correlation functions are given by: 



(n(i)n a {i)) 



1 



-{G 1 X 1 +aG 2 X 2 [l + X 1 



M*K(«)> 



(D(i)D a (i)) = 



+2b(X 1 - 2X 3 )} + (Gi - 2G 3 ) 
[a 2 X 2 + 6(1 + X l ){X l - 2X 3 ) + 
+b 2 (X 1 -2X 3 ) 2 }} , (27) 

{G2X2 + G 2 X 2 [b + (a 2 + 

(e pH i ) 

+b + b 2 )(X 1 -2X 3 )]+a(G 1 -2G 3 ) 

[(l + b)Xl + +(X 1 -2X 3 ) 

[l + b(X 1 -2X 3 )}}} , (28) 

G * X * (29) 



The internal energy per site is given by: 

E = UD - J (n 3 {i)n^(i)} - 2hm . (30) 

Specific heat G, charge \c and spin Xs susceptibilities are 
given by: 

dE dn dm 

c= dr> Xc = T d; ' Xs= m- (31) 



B. Half- filling case 

When n = 1, equation (|2"2"|) gives Gi = 1. Recalling 
(IB161) : 



(32) 



in agreement with the particle-hole requirement fi(2 — 
n) = U — fJ,(n). 

Then, the self-consistent equations (|B24[) - (IB26|) take 
the form: 

Fi = (l-X 1 )[l + aG 2 X 2 + 

+6(l-2G 3 )(Xi -2JC 3 ) =0, (33) 

F 2 = G 2 -X 2 + (l-2G 3 )X 2 {2a-b + 
- [b + {a - bf] (Xt - 2X 3 ) + 
+G 2 {aX 2 (a -b-l) + (Xi- 2X 3 ) ■ 
■\2b-a + b{b-a){Xx -2X 3 )]} =0 , (34) 

F 3 = G 3 - X 3 [l + aG 2 X 2 + 6(1 - 2G 3 ) ■ 

■(X 1 -2X 3 )]=0. (35) 



The first equation has the solutions: 
X 1 = 1 

■y {l+aG 2 X 2 ) , 0Y 
Xl 6(1-2G 3 ) + 2X 3 



(36) 



For the second solution, substituting in (1351) we obtain 
G 3 = which is clearly not satisfied. Therefore, for n = 1 
we have X\ = 1 , and we remain with two parameters X 2 , 
X 3 to be determined by means of the equations: 

F 2 = G 2 -X 2 + (l-2G 3 )X 2 {2a-b + 
- [6+(a-6) 2 ] (1-2X3) + 
+G 2 {aX 2 (a - b - 1) + (1 - 2X 3 )- 
• [26 - a + 6(6 - a)(l - 2X 3 )}} = , 

F 3 = G 3 - X 3 [1 + aG 2 X 2 + 6(1 - 2G 3 ) • 
•(I-2X3)] = 0. 

C. Density of states 

By noting that the cross GFs (R (i, t)ip ^ (i, t')] ) 
vanish, the electronic density of states (DOS) for spin a 
is given by 



Im 



2z+l 

E 



g^h + g^h 



Therefore we have: 

/-(-00 2z+l - 

^x CT (w) = 
-°° n=l 



(£,") , (»7,n) 
Per, 11 + AVll 



= 1 



(37) 



(38) 



D. Isomorphism with a two- level Ising model 



Let us consider the following transformation: 
nt(<) = \ [1 + , n 4 (i) = i [1 + &(<)] . (39) 

From the properties of n^(i) and 714. (i), it is clear that 
Si(i) and are two spin variables with S 1,2(1) = ±1. 
According to this transformation the Hamiltonian @ 
reads as: 



H = J/s^SitySZii) + S2(i)S2(i)) + 

i 

+ Jj. + J± ^2 Si(i)S$(i)+ (40) 

i i 
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and describes two interacting spin chains in the presence 
of an external magnetic field. Here J// = — J/4 is the 
NN intra-chain interaction, J± = {7/4 and J" = J/2 are 
NN and next-NN inter-chain interactions, respectively. 
Importantly, in both chains the spins are coupled with 
an effective chain-dependent magnetic field h n — {7/4 — 
fjt/2 + (-l) n h/2 (n = 1,2). In particular, at half-filling 
we have /j, = U/2 and an alternating magnetic field h n = 
(-l) n h/2. 



E. Summary 

Summarizing, in this Section we have shown that the 
atomic Hubbard model with magnetic inter-site interac- 
tions is in principle exactly solvable for any dimension of 
the system by using the Green's function formalism and 
the equation of motion method. The central point is the 
fact that there exists a closed set of eigenoperators of the 
Hamiltonian ((2} [see Eqs. (j9|)] allowing to determine ex- 
act expressions of the Green's and correlation functions. 
These expressions are determined in terms of few local 
correlators and . For the one-dimensional case, 
by using algebraic relations of the relevant operators and 
symmetry properties, these correlators can be expressed 
in terms of the chemical potential /i and of the basic 
parameters X\, X2, and X3. These four quantities are 
determined as functions of the external parameters n, U , 
J, h, T by solving the system of equations (|B24l) - ([35)l and 
(|2"2l . Once these quantities are known, all properties of 
the model can be computed. In the next Section we shall 
present a comprehensive study of the zero-temperature 
phase diagram and report the distribution of the state 
density by identifying the different phases in the plane 
(U,h). We shall also study the behavior of several quan- 
tities such as specific heat, charge and spin susceptibil- 
ities, and density of states, in the limit of vanishing T. 
A detailed analysis on the temperature dependence of all 
the aforementioned quantities is reported in Sec lIVl 



III. ZERO TEMPERATURE LIMIT 

In this Section we present several results obtained in 
the limit of zero temperature. Hereafter we will put 
I J| = 1 in order to fix the energy scale and analyze both 
ferro (J = 1) and antiferro (J = — 1) inter-site magnetic 
couplings. 



due to the invariance under h — » —h, n-f — > nj, trans- 
formation, we restrict our analysis to h ^ . One can 
distinguish two different cases: J = 1, J = — 1. 



1. Ferromagnetic inter-site coupling 

The phase diagram at zero temperature for positive 
values of J is shown in Fig. [TJ where we consider the 
(U,h) plane. Two different phases can be identified. 




Figure 1: 
T = 0. 



(Color online) 2D phase diagrams for J = 1 at 



NM-phase. This phase, which we call "non-magnetic" 
(NM) phase, is observed for U <U C where U c = — 2( J + 
h) and for ^ n ^ 2. It is characterized by the absence 
of magnetization and spin-spin correlations; there is no 
dependence on the external magnetic field h. The attrac- 
tive local potential U prevails on both the magnetic field 
and the ferromagnetic coupling; there are no sites singly 
occupied and the double occupancy D is an increasing 
function of the filling. The chemical potential depends 
only on U; the electrons are not correlated and the two- 
site charge and double occupancy correlation functions 
are the squares of the corresponding one-site ones: 



fi = U/2 
(£>(<)) = n/2 
(m(i)) = 



(n(i)n a (i)) = 

(n 3 (i)nf(z)) 

(D(i)D«(i)) 





(n/2) 2 



(41) 

Recalling (|30| . it is immediate to see that the internal 
energy per site has the value Enm = nU/2. A typi- 
cal configuration occurring in this phase is a mixture of 
doubly occupied and empty sites and is shown in Fig. EI 



ti*l 



A. Phase diagram 

In the limit of zero temperature we expect to find 
charge and spin long-range orders. To characterize all 
the possible configurations, we have solved numerically 
the set of equations (IB24I) - (|B26|) . (|2"2"j) and studied the 
trends of each relevant correlation function depending on 
the values of n, J, h and U. For the sake of simplicity, 



Figure 2: One of the possible spin and charge configurations 
for the NM-phase at T = 0, J = 1. t an d 4- represent the two 
possible spin states. 

F-phase. This phase, called "ferromagnetic" (F) phase, 
is observed for U > U c and for ^ n ^ 2. It is character- 
ized by a dominant ferromagnetic order and finite values 
of magnetization and spin-spin correlation functions. In 



the region ^ n ^ 1 the double occupancy is zero; all 
the spins are aligned and the magnetization reaches its 
saturation value (m(i)) = n/2; the spin-spin correlation 
function is positive and increases with n. The chemical 
potential does not depend on the filling and is a function 
of J and h: 



|U = — J — h (n(i)n a (i)} = n 
(D(i)) = Q (n 3 (iK(i)) = n 
(m(i)) = n/2 (D(i)D a (i)) = 



,0 < n < 1 



(42) 



At n — 1 the chemical potential exhibits a discontinuity 
and jumps at fj, = U/2 . In the region of filling n > 1 the 
chemical potential takes the value a — J+h+U, in agree- 
ment with the scaling law fi(2 — n) = U — u{n). The dou- 
ble occupancy increases linearly with n; correspondingly, 
the magnetization and the spin-spin correlation function 
{n3(i)n^(i)) decrease: 

H = J + h + U (n{i)n a (i)} = 3n - 2 
(D(i))=n-1 (n 3 (i)n%(i)) = 2 - n , 1 < n < 2 . 
(m(i)) = 1 - n/2 (D(i)D a (i)} = n - 1 

(43) 

The internal energy has the value: 

f-n(J + ^) ,0<n<l 

1 U(n-l)-(2-n)(J + h) , Un<2 v 7 

A typical configuration occurring in this phase is shown 
in Fig. |3] for the cases n < 1 and n > 1. We see that: 

f-n(j + ft+§) ,0<n<l 
^-• B ^ = j- (2 -n)(J + /i+f) ,1<„<2' 

(45) 

therefore, regardless of the specific value of filling n, there 
is a critical value of the local potential U c = — 2( J + h) 
which separates the two phases. 




Figure 4: (Color online) Signature of NM-F phase transi- 
tion in magnetization (left) and spin-spin correlation func- 
tion (right) plotted as functions of the magnetic field h, for 
different values of U at n = 1 and J = 1. 



2. Antiferromagnetic inter-site coupling 

The phase diagram at zero temperature for negative 
values of J is shown in Fig. [5j where we consider the 2D 
(U,h) and the 3D (U,h,n) cases. 



J=-1.0 ; 0.0<n<0.5 
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Figure 3: Some of the possible spin and charge configurations 
for the F-phase at T = 0, J = 1. 

Signatures of the NM-F phase transition occurring in 
magnetization and spin-spin correlation function are re- 
ported in Fig. [4] 

It is immediate to see that for U < —2, there is a 
phase transition from the NM phase to the F phase at 
h = h c = —J — U/2. According to (jlTj) . both mag- 
netization and spin-spin correlation functions are zero 
for h < h c . Vice-versa, when h > h Cl (m(i)) and 
(713(1)713 (z)) assume finite values according to P2"j) -(f4*3" |) . 
For U ^ —2 instead, there is no phase transition: the sys- 
tem is always in the F phase. In particular, for zero field 
there is a magnetic order: the spin-spin correlation func- 
tion is equal to 71 while the magnetization jumps from 
—n/2 to n/2 crossing h = 0. 



Figure 5: (Color online) 2D and 3D phase diagrams for J = 
-1.0 at T = 0. 

In the regions ^ n ^ 0.5 and 1.5 ^ n ^ 2 the phase 
diagram does not depend on n (as shown in Figs. [5ji 
and [3>) and is characterized by three different phases 
which join at the tri-critical point P\={U = —2, h = 1}. 
In the region 0.5 ^ n ^ 1.5 the phase diagram depends 
on 11 and is characterized by four different phases and 
two tricritical points: P\ and P2, as shown in Figs. Eh 
and [5]i. It is worth noting that while Pi is fixed, the 
position of P2 in the U-h plane changes as a function of 
n, as reported in Fig. [7l We can also note that, changing 
the filling, Pi and P2 never coincide and remain always 
well separated from each other. 

NM-phase. Also for negative values of J a "non- 
magnetic" (NM) phase is observed, characterized by the 
absence of magnetization and spin-spin correlations. All 
the features of this phase have already been reported in 
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Eqs.gJ). 

AF-phase. For U > —2 J and for low values of the mag- 
netic field h an antiferromagnetic (AF) order is observed 
in the entire region of filling. This phase is characterized 
by the absence of magnetization and negative values of 
spin-spin correlation function (713(2)713 (i)): 



(46) 

For ^ 71 ^ 1, as the filling increases, the electrons 
start to singly occupy the available empty sites aligning 
their spins opposite with respect to their nearest neigh- 
bors; correspondingly, the magnetization and the double 
occupancy are zero. The half-filling configuration corre- 
sponds to a perfect Neel state. For n > 1, by increasing 
the filling, some of the sites become doubly occupied, 
while the magnetization remains zero. A typical config- 
uration which takes place in this phase is shown in Fig. 
[5] for the cases n = 0.5, n = 1 and n = 1.5. 
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Figure 6: Some of the possible spin and charge configurations 
for the AF-phase at T = 0, J = — 1. 



Fl-phase. When the external magnetic field is strong 
enough to dominate with respect to J, a ferromagnetic 
behavior (Fl) is observed despite of the presence of an 
anti-ferromagnetic coupling. In this phase, the effect of 
the magnetic field induces the spins to be all aligned and 
the magnetization reaches its saturation value. Due to 
the presence of an anti-ferromagnetic coupling, contrarily 
to what has been said for the F-phase, NN spin-spin and 
charge-charge correlation functions remain zero until n = 
0.5. At n = 0.5 a charge order state is observed with a 
checkerboard structure. 
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Figure 7: (Color online) Position of P2 3-critical point in the 
h-U plane (left) and in the 3D phase diagram (right). 

On the contrary, as reported in Fig. HI for 0.5 < 
71 < 1.5 the spin-spin correlation function becomes finite; 
(713(1)713(1)) increases with n up to half-filling, then de- 
creases and vanishes for n 1.5. For 0.5 < n ^ 1.0, D — 
and the sites are singly occupied by electrons with par- 
allel spins. Exactly at half-filling, all the sites are singly 
occupied with polarized spins, therefore the magnetiza- 
tion assumes its maximum value. For 1.0 ^ n < 1.5, 
the double occupancy becomes finite but the NN cor- 
relation function (D(i)D a (i)) is always zero; exactly at 
71 = 1.5 another checkerboard structure with a pattern 
composed of alternating singly and doubly occupied sites 
is observed. For n ^ 1.5, {D(i)) and (D(i)D a (i)) increase 
up to the maximum value, while the magnetization goes 
to zero as n — > 2. 
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Figure 8: Some of the possible spin and charge configurations 
for the Fl-phase at T = 0, J = -1. 
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(47) 



F2-phase. As shown in Fig. [5l in the range in which /i 
and J are comparable, for J7 < and for 0.5 < n < 1.5 
the model exhibits an anomalous ferromagnetic phase 
(F2) induced by the competition between the magnetic 
field and the antiferromagnetic inter-site coupling. This 
phase is characterized by the absence of NN spin-spin 
correlations and the presence of a constant magnetiza- 
tion. 
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/i = U/2 (n(i)n a {i)} =2n-l 

(£>(*)) = n/2 - 1/4 (ns(<K(<)>=0 . 

(m(i)) = 1/4 (D(i)D a (i)) = ^ 
S F2 = C/(|-|)-| 

As shown in Fig. [HI in this phase both (D(i)) and 
(n(i)n a (i)) are finite while (713(1)713 (i)) remains zero. 
At n = 0.5, doubly occupied sites appear in between 
two singly occupied ones; by increasing n, the number 
of doubly occupied sites increases, while the number of 
singly occupied sites remains constant. This explains 
why the magnetization does not change. The n = 1 
configuration is characterized by a particular pattern in 
which one or more clusters, composed of the following 
pattern: (cr,t^,cr,0), are repeated periodically along the 
chain. Some signatures of NM-F1-F2 and AF-F1-F2 
phase transitions occurring in magnetization and spin- 
spin correlation function are reported in Fig. 1101 



[n^0.25) 



[n^0.5] 







u 






t \n\ t 




t u 








U 








A v 






































[n=a.5J 


A 


u 




n 


t |u| t 




t u 




N 


A 


U 




A v 


A 


A v 


t u 




































[n=a.75J 


u 


u 


n 


'■ 


u|u|u 




tin 


u 


A 




u 


\ 


A 


\ 


A v 


u t 



Figure 9: Some of the possible spin and charge configurations 
for the F2-phase at T = 0, J = -1. 



is not observed: the magnetization jumps from zero (AF- 
phase) to n/2 (Fl-phase) while the spin-spin correlation 
function goes from negative (AF-phase) to positive (Fl- 
phase) values. 



B. Charge and spin excitations 

The presence of different phases and long-range orders 
at T — >• can also be predicted by looking at charge and 
spin excitations in system's response functions such as 
charge and spin susceptibilities. As shown in the previous 
Subsection, in each phase at T — all single particle 
properties and correlation functions depend only on n but 
not on the model parameters (U, J, h). Therefore charge 
and spin susceptibilities are expected to have a constant 
value in all the phases with divergencies or discontinuities 
localized at the phase boundaries where small variations 
of the model parameters can imply the transition from 
a charge/spin ordering to another. In particular, while 
Xs{T = 0) = in each phase, for the charge susceptibility 
we have instead: 



Phase 
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n = 0, or, n = 1 


Fl 


n(n - l)(2n - 1) 
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2rc(2-n) -3/2 


0.5 < n sC 1 
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Figure 10: (Color online) Signatures of phase transition in 
magnetization (left) and spin-spin correlation function (right) 
plotted as functions of the magnetic field h, for different values 
of U at n — 1 and J = — 1. 

As shown in Fig. [TQl for U ^ —2, changing h one 
crosses three different phases: NM, Fl and F2. In this 
range (solid curves) the magnetization assumes three dif- 
ferent values: in the NM-phase, 1/4 in the F2-phase 
and n/2 in the Fl-phase. The spin-spin correlation func- 
tion jumps from zero (in NM and F2 phases) to one (in 
the Fl-phase). Instead, at — 2 < U < (dashed curves), 
NM-phase is replaced by AF-phase and the spin-spin cor- 
relation function becomes negative when \h\ $J 1.5. Fi- 
nally, for positive values of U (solid curves) the F2-phase 



As reported in Fig. [IT] and Fig. [T^J the zero tem- 
perature behavior of charge/spin susceptibilities shows a 
discontinuity /divergence crossing the phases character- 
ized by different charge/spin orderings, reproducing the 
boundaries of the phase diagram of the model for both 
J = 1 and J = — 1. 

Importantly, it is worth noting that the phase diagram 
of the model can be reconstructed by joining together all 
the information coming from both charge and spin sus- 
ceptibilities since different phases can exhibit the same 
charge ordering but different spin orderings. In this re- 
gard, we note for example that at n = 1 the charge sus- 
ceptibility is not able to distinguish between AF and F2 
phases since both of them at n — 1 are characterized by 
the same charge ordering with all sites singly occupied. 
Similarly, spin susceptibility is not able to distinguish 
between NM and AF phases since both of them are char- 
acterized by the absence of magnetization. 

The same analysis can also be done in terms of the 
specific heat C(T) whose behavior in the T — » limit, as 
we shall see in detail in the next Section, is characterized 
by low-temperature peaks in the proximity of the phase 
boundaries. This analysis does not allow to distinguish 
between charge and spin excitations; however, as shown 
in FigHSl it provides a good estimate of all the phase 
boundaries of the model. 
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Figure 11: (Color online) Contour-plot of the charge suscep- 
tibility at T = 0.001 as a function of the external magnetic 
field h and the local Coulomb potential U for: (a) J — 1, 
n = 0.5; (6) J = 1, n = 1; (c) J = —1, n = 0.75; (d) J = -1, 
n = 1. 





Figure 12: (Color online) 3D plot of spin susceptibility in the 
limit of zero temperature (T = 0.001) at n = 1 as a function 
of the external magnetic field h and local Coulomb potential 
U for J = 1 (left) and J = -1 (right). 



C. Density of states 

To complete the zero temperature analysis, we report 
in this Subsection the results obtained for the density 
of states in the limit of zero temperature. As shown in 
Eg. (1571) . the density of states is expressed as a superpo- 
sition of several delta functions, each of them centered 
at the energy levels and weighted by the spectral 
functions p( a < n > . In the limit of zero temperature, most 
of the weights vanish and only few energies, correspond- 
ing to the ground state and first excited states, give a 
contribution. We report in Fig. Q3] and Fig. [TS] the den- 
sity of states calculated for each phase that appears in 
the J = 1 and J = — 1 phase diagrams at half- filling and 
T = 0.001. 

For ferromagnetic inter-site coupling, as shown in Fig. 
[31 the NM-phase at n = 1 is characterized by the presence 



Figure 13: (Color online) Contour-plot of the specific heat at 
low temperature at half filling for J = 1 (left) and J = — 1 
(right) as a function of the external magnetic field h and local 
Coulomb potential U. 
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Figure 14: (Color online) Density of states in the limit of zero 
temperature for the NM and F phases ( J = 1) . We report 
the total density of states (c(i) = £(i) + r](i)) contributions 
for both spin up (a — +1) and spin down [a = —1). The 
contributions due only to and r)(i) fields are also reported 
in the insets. 



of either doubly occupied or empty sites. Accordingly, 
the low-laying excited states contain singly occupied sites 
with spins aligned along h. Therefore, as shown in Fig. 
Q3J the density of states has a peak above the Fermi level 
induced by £f(i). The F-phase at n = 1 is composed of 
singly occupied sites with spins pointing towards h. On 
the contrary, the low-laying excited states admit a num- 
ber of doubly occupied sites, resulting in a peak above 
the chemical potential induced by rj^, as shown in Fig. 

M 

The same considerations can also be done in the case 
of negative J (see Fig. [15]) so that the nature of the peaks 
around the Fermi level can be easily predicted by observ- 
ing the zero-temperature configuration of each phase. It 
is worth noting that, only in the F2-phase, both the con- 
tributions of £ and rj are allowed thanks to the fact that 
any F2-configuration comprises either polarized spins or 
double occupancies (Fig. [§]). 

Our analysis of the density of states also allows to get 
some insight about the gap A that separates the ground 
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Figure 15: (Color online) Density of states in the limit of zero 
temperature for the AF, Fl and F2 phases (J = — 1). We 
report the contributions for both spin up (a = +1) and spin 
down [a — — 1). In the insets the contributions due to 
and rj(i) fields are reported. 



state from the first excited one. It can be seen that, as 
long as the system remains in the same phase, a finite gap 
exists according to the results obtained for the specific 
heat in the limit of zero temperature (Fig [13]). On the 
contrary, moving from a phase to another, A closes as one 
approaches the phase boundaries and low-temperature 
excitations appear in the specific heat signalling a level- 
crossing between the ground and the first excited state. 



D. Comparison with variational approach 

In closing this Section, to further emphasize the crucial 
role played by our exact solution as a guideline for test- 
ing analytical and numerical approximations, we com- 
pare our results with those obtained in Ref|H, where the 
same model, but in absence of magnetic field, has been 
studied within a variational and a mean-field approach. 
First of all, the absence of magnetic field introduces a 
symmetry between the ferro (J > 0) and the antiferro 
(J < 0) case. Owing to this reason, only the case J > 
has been studied in Refl48l. where the ground state phase 
diagram exhibits only the two phases NM and F (or AF 
for J < 0). We note that, as a consequence of the mean- 
field approach the ground state energy in the F phase 
(called EF Refl48T) reported in RefJH does not agree with 
ours, given in Eq. (Pffl . In particular EF of Refl48l de- 



pends on the square of the filling and this implies that 
the transition line between the NM and F phases depend 
on n in strong contrast with our exact calculations re- 
ported in Eq05j Another important difference is that the 
thermodynamic instability of the NM phase, predicted in 
RefHa. is not observed in our exact solution. 



IV. FINITE TEMPERATURE RESULTS 

In this Section we present the study of the model ([2]) 
at finite temperature. As already done in the previ- 
ous Section, we put |J| = 1 and analyze both J = 1 
and J = — 1 cases. We report the temperature depen- 
dence of single-particle (chemical potential, double oc- 
cupancy, magnetization) and thermodynamic (entropy, 
specific heat, charge and spin susceptibility) quantities 
for different values of n, U, h and J. Because of the 
one-dimensionality of the model, there is no long-range 
order at finite temperature. However, in this Section we 
maintain the reference to the T = phases, described in 
Section III, in the sense that we label each set of param- 
eter (n, U, h, J) according to the corresponding phase 
observed at zero temperature. 

Chemical potential. The behavior of the chemical po- 
tential /i as a function of the particle density n is reported 
in Fig. [16] for different values of T, J, U, and h. In the 
magnetic phases, an interesting feature is the presence 
of crossing points in the chemical potential curves, when 
plotted versus the filling for different temperatures. More 
precisely, one observes crossing points at the commen- 
surate fillings which also are turning points where the 
derivative dfj,/dT changes sign. 
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Figure 16: (Color online) Chemical potential plotted as a 
function of the filling for different temperatures. Different 
values of h and U are related to NM (top, left), F (top, right), 
Fl (bottom, right) and F2 (bottom, left) phases. 

Double occupancy. We report in Fig. [17] the tempera- 
ture dependence of the double occupancy at n — 1 and 
different values of U and h, for both J = 1 and J = — 1 



12 



cases. 

In the NM phase (Fig. [T7]a, U = -4 and U = -3), the 
attractive local potential U prevails and at zero tempera- 
ture D reaches its maximum value D = n/2; as the tem- 
perature increases, the thermal excitations tend to break 
the doublons (the pairs of electrons with anti-parallel 
spins residing on the same site); several configurations 
with singly occupied sites become available resulting in 
a depression of the double occupancy. D appears always 
a decreasing function of the temperature, going to n 2 /4 
in the limit T — > oo. 

In the F phase (Fig. Ek, U = -2.8 and U = -2), 
the magnetic interactions prevail on the attractive local 
potential U and D vanishes at T = 0. For U > U c two 
behaviors can be distinguished. When the temperature 
is raised above zero, D increases: because of the presence 
of an attractive on-site U potential, the thermal fluctu- 
ations favor the formation of doublons, and the relative 
contribution of the excited states with large average D 
increases. By further increasing T, the double occupancy 
D reaches a maximum, then decreases. The temperature 
where D exhibits a maximum decreases with U and goes 
to zero for U — > U c . When U is very close to U c (e.g. 
U = —2.9) at the temperatures of the order of the gap 
between the first excited and the ground state, the two 
become quasi-degenerate, and D acquires a contribution 
from the latter. In the limit of U — > U c , close to T = 
the double occupancy exhibits a jump from zero to n/2. 

Similar considerations can be done for the J = — 1 case. 
As reported in Fig. [Tfb . in AF (ft, = 0.1 and h = 1.4) and 
Fl [h = 2.6 and h — 4.0) phases double occupancy is sup- 
pressed by the presence of magnetic orders at T = 0. By 
increasing the temperature, D monotonically increases. 
On the contrary, in the F2-phase (h = 2.0 and 2.4), de- 
spite the presence of a ferromagnetic order, one has a 
finite double occupancy at T = [D(0) = (2n - l)/4], 
resulting in a non-monotonic T-dependence. In all phases 
D tends to n 2 /4 in the limit of T — > oo . 
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Figure 17: (Color online) The double occupancy as a function 
of temperature for J = 1, h = 0.5, n = 1, —4 ^ U ^ —2 (left 
panel) and J = -1, U = -1.0, n = 1, 0.0 s£ h < 4.0 (right 
panel) . 

Magnetization. We report in Fig. [18] the temperature 
dependence of the magnetization at different values of U 
and h, for both J = 1 and J = — 1 cases at n = 1. In the 
NM phase (Fig. [Tgk, U = -4 and U = -3) the ground 



state is composed of double occupied sites only (see Fig. 
[2]). By increasing T, excitations destroy doublons and a 
finite magnetization appears. The increase of the magne- 
tization with temperature is contrasted by thermal fluc- 
tuations which eventually destroy m at T — > oo. In 
the F-phase (Fig. [Tgja,, U = -2 and U = -2.8) the 
spins in the ground state are fully polarized. Therefore 
m(T = 0) = n/2. Increasing the temperature, doubly 
occupied sites progressively emerge, lowering the magne- 
tization that also exhibits a broad maximum at interme- 
diate temperatures associated with higher energy excited 
states. 

Apart from the NM-phase, in the case of antiferro- 
magnetic coupling, three different behaviors of magne- 
tization can be distinguished which correspond to three 
different phases. As reported in Fig. [SJ in the Fl-phase 
at T = all the spins are aligned in the field direction 
and m(T = 0) assumes its saturation value. Thermal 
fluctuations destroy the magnetization which decreases 
monotonically upon the increase of temperature. In the 
AF phase the average magnetization per site is zero at 
T = because of the antiferromagnetic ordering. At fi- 
nite T, the magnetization acquires contributions from the 
excited states which are magnetic only at finite h. There- 
fore, for h ^ the magnetization exhibits a maximum at 
an intermediate temperature T* and goes gradually to 
zero as T — > oo. As shown in Fig. |H1 in the F2-phase 
(h = 1.6 and h — 2.4) the ground state is characterized 
by the coexistence of singly and doubly occupied sites 
under the constraint that m = 1/4 independently on the 
value of n € [0.5, 1.5]. Upon the increase of T, thermal 
excitations break some of the doublons and misalign the 
spins. Hence there is a high competition among multiple 
energy scales: T, h, J, U and the magnetization may 
exhibit one or several peaks, each of them corresponding 
to a distinct energy scale. 
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Figure 18: (Color online) The magnetization as a function of 
temperature for J = 1, h = 0.5, n = 1, —4 ^ U ^ —2 (left 
panel) and J = -1, U = -1.0, n = 1, 0.0 s£ h < 4.0 (right 
panel) . 

Internal energy. In the NM phase, as shown in Fig. 
119b .. the internal energy monotonically decreases by low- 
ering T and tends to the value Enm = nU/2 as T — > 0. 
In Fig. [T9b . the temperature dependence of the internal 
energy E is shown for the F-phase, by considering various 
values of U and n = 1. E decreases with T and tends to 
the value Ep — —n(J + h) as T — > 0. When U is close to 
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U c , E exhibits a discontinuity at a certain temperature, 
which goes to zero as U tends to U c . 

For the case J = — 1, interesting features for the inter- 
nal energy are observed (see Fig. IT5b ) only for the AF- 
phase in proximity of the transition to the NM-phase. 
For all the other phases, the energy is found to be a con- 
tinuous and increasing function of the temperature. The 
significant variations observed in the E(T) dependence 
lead to the presence of very low temperature features in 
the specific heat, as we shall discuss below. 
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Figure 19: (Color online) Internal energy as a function of 
temperature for NM-phase (panel a) and F-phase (panel b) 
at J = 1 and AF- phase (panel c) at J = — 1. 

Charge susceptibility. We report in Fig. [20] the tem- 
perature dependence of the charge susceptibility \c at 
different values of n, U and h for the J = 1 case. In the 
NM phase (Fig. l20k , U = -4 and U = -3), Xc takes the 
value Xc = n{2 — n) at T = 0, independently on the value 
of U, and decreases by increasing T . For a given tem- 
perature, Xc decreases by increasing U and by decreasing 
the filling. In the F phase at n = 1 (Fig. EOK U = -2.8 
and U = —2) Xc vanishes at T — 0, increases with T up 
to a maximum and tends to a finite value in the limit 
T — > oo. When U approaches U c , the maximum becomes 
more pronounced and tends to one as T — > 0. As shown 
in Fig. [20b . in the F phase Xc has a different behavior 
away from half filling. The charge susceptibility increases 
by lowering T and tends to diverge as T —y 0; when U 
approaches U c , the tendency to diverge becomes much 
stronger. In the limit of high temperatures, for both the 
NM and F phases, the charge susceptibility tends to a 
constant value which does not depend on U but only on 
n according to the law: limT->oo Xc = n(2 — n)/2. 




Figure 20: (Color online) The charge susceptibility as a 
function of temperature at J = 1, h = 0.5 for: n = 1, 
-4 U < -2 (left panel) and U = -2 and 0.4 sC n sC 1 
(right panel). 




Figure 21: (Color online) The charge susceptibility as a func- 
tion of the filling and different temperatures for AF-phase (left 
panel) and Fl (central panel) and F2 phases (right panel). 

In the case of anti- ferromagnetic coupling (J = — 1), 
the AF phase exhibits the same temperature and filling 
dependence already described for the F phase. As re- 
ported in Fig. l2"Th . for any n / 1 the charge susceptibil- 
ity increases as the temperature decreases and diverges 
for T — > 0. At half-filling Xc has a peak which gradually 
moves towards T — as the system approaches the phase 
boundary. 

On the contrary, as reported in Figl2Tb and Figl2Th. in 
Fl and F2 phases there are no divergencies in the charge 
susceptibility at any value of filling and temperature. In 
the Fl-phase fFig[2"Tb) Xc(T — 0) is zero at any value of 
the filling which corresponds to the presence of ordered 
phases at T = (n = 0,1/2,1,3/2,2). In the F2-phase 
fFig ]21b ) Xc has a maximum at n = 1, then decreases as 
one approaches to n = 1/2 and n = 3/2. These values 
of filling are situated at the phase boundary between the 
Fl and F2 phases, therefore Xc{T = 0) = in agreement 
with the above discussion. It is worth noting that, sim- 
ilarly to the J = 1 case, limr->oo Xc = n(2 — n)/2 in all 
the phases. 

Spin susceptibility. In Fig. 122b . the temperature de- 
pendence of the spin susceptibility Xs is shown for the 
NM phase, by considering J=l, h = 0.5, n = 0.8 and 
various values of U. Xs vanishes at T = 0, increases by 
increasing T up to a maximum temperature T max , then 
decreases. This behavior is observed for all values of the 
filling. For a given temperature, Xs decreases by decreas- 
ing U and by decreasing the filling. The temperature 
Tmax decreases by increasing U. In Fig. [22b . the tem- 
perature dependence of the spin susceptibility is shown 
for the F-phase, by considering various values of n and 
U = —2.8. Xs vanishes at T = Q, increases by increas- 
ing T up to a maximum temperature T s , then decreases. 
When U approaches U c , the temperature T s moves to- 
ward zero and Xs tends to diverge. 

In the case of anti- ferromagnetic coupling (J = — 1), 
Xs{T) in the AF-phase (see Fig. HSi, U = -1.9 and 
U = —1.0) exhibits the same features already analyzed 
for NM and F phases with a maximum that moves to 
T = as the phase boundary is approached. The behav- 
ior of Xs{T) in the Fl and F2 phases is shown in Fig. 
123b . Low-temperature peaks, whose heights increase by 
approaching the transition line, are observed in the Fl 
phase and are due to the contribution from the low-laying 
magnetic excited states. 

Concerning the filling dependence, we can note that 
in all the phases Xs(T = 0) =0 in the whole range of 
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Figure 22: (Color online) The spin susceptibility \s as a func- 
tion of T for J — 1, h = 0.5; left panel: n = 0.8 and various 
values of U ; right panel U — —2.8 and various values of n. 
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Figure 23: (Color online) The spin susceptibility \s as a func- 
tion of T for J = — 1, n = 1.0 and various values of U for 
h — 3.0 (left panel) and h = 0.5 (right panel). 



filling. By increasing the temperature the spin suscep- 
tibility also increases with different temperature depen- 
dencies according to the particular choice of the exter- 
nal parameters. In the AF phase, as shown in Figl24k. 
Xs 7^ at any n ^ 1 at low-temperatures, with a max- 
imum at n — 1/2, 3/2 moving to n = 1 with increas- 
ing T. In the Fl-phase instead, as shown in Fig[2"4"b. 
Xs {T 7^ 0) exhibits a very narrow low-temperature peak, 
centered at n = 1, whose dispersion increases with in- 
creasing temperature. This feature also clearly appears 
in the crossover from Fl to F2 phase. In fact, as reported 
in Fig|24b. Xs(T) remains finite in a wide range of tem- 
peratures for 1/2 ^ n ^ 3/2 (F2-phase) while rapidly 
decreases with decreasing temperature in the Fl-phase 
for sC n < 1/2 and 3/2 < n ^ 2. 

It is worth noting that, independently on the sign of 
the spin exchange J, in the limit of high temperatures the 
spin susceptibility follows the Curie law with a coefficient 
which does not depend on U but only on fi\ lim^—^oo Xs — 
n(2 — n)/AT. Therefore we have limi^oo Xc/Xs = 2T, in 
which the factor 2 is due to the spin multiplicity. 

Specific heat. In the NM phase the specific heat ex- 
hibits a peak at high temperature (see Fig. [23k). The 
position of this peak, T2, and its intensity, hi, decrease 
by increasing U . T2 decreases with almost a linear law. 
At finite temperature, the possible excitations are transi- 
tions to configurations where some sites are singly occu- 
pied, with spins aligned in the direction of the magnetic 

(ri) 

field. The field ipy is responsible for these transitions 
and this explains why T2 is linear in U . It is worth noting 




Figure 24: (Color online) The spin susceptibility as a function 
of the filling and different temperatures for AF-phase (panel 
a), Fl-phase (panel b) and Fl (0 < n < 1/2, 3/2 < n < 2) - 
F2 (1/2 < n < 3/2) phases (panel c). 



that the temperature T2, at which the specific heat has 
a maximum, is close T m /2, where T m is the position of 
the peak in the magnetization. 

The specific heat in the F-phase is shown in Fig. [25b . 
In general, in the F-phase C exhibits lower (Ti) and 
higher (T2) temperature peaks, associated to transitions 
involving the ground state and the low-laying excited 
states. Approaching the boundaries, the ground and the 
first excited states become quasi degenerate. Hence a 
temperature of the order of the gap between these two 
states will be sufficient to induce thermal excitations, re- 
sulting in the presence of a low-temperature peak that 
moves towards T = 0. On the contrary, position and in- 
tensities of high-temperature excitations, responsible for 
a second lower and broadened peak, remain almost un- 
changed. In fact, as shown in Fig. [26b . the intensity 
/12 of the higher temperature peak remains almost con- 
stant when U varies. Contrarily, h\ rapidly increases 
as U approaches U c and diverges in the limit U — > U c . 
This occurs since the degeneracy of the first excited state 
(belonging to the NM-phase) is infinite in the thermody- 
namic limit with respect to the degeneracy of the ground 
state (F-phase). 





Figure 25: (Color online) Specific heat as a function of tem- 
perature for J = 1, h = 0.5, n = 0.8, -4 < U < -3 (left 
panel) and —3 ^ U ^ —2 (right panel). 

Most of the features observed in the F-phase persist for 
J = — 1 in the AF-phase. As shown in Fig. [27k . the spe- 
cific heat exhibits a low-temperature peak whose height 
(dispersion) increases (decreases) as one approaches the 
phase boundary. As reported in Fig. [51 the AF- 
configuration at n — 1 is composed exclusively of singly- 
occupied sites. Therefore, the low-laying features cor- 
respond to transitions from the ground states to some 
excited states with finite D or m ^ as confirmed by 
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Figure 26: (Color online) Positions (Ti in F-phase and T2 
in NM-phase) and intensities (/ii in F-phase and /12 in NM- 
phase) where the specific heat has maximas as functions of U 
for J = 1, h = 0.5, and different values of the filling. 
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Figure 27: (Color online) Specific heat as a function of tem- 
perature at J = —1 for AF (left panel), Fl (central panel) 
and F2 (right panel) phases at n = 1.0. 



the analysis of the temperature dependence of double 
occupancy (FigfTT]) and magnetization (FigHTf. On the 
contrary, Fl and F2 phases are characterized by both low 
and high temperatures peaks whose heights remain quite 
constant when the phase boundaries are approached. 
While the position of the former remains roughly con- 
stant, the latter moves towards T — 0. 
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Figure 28: (Color online) Low-temperature features of the 
specific heat at different transition lines for J = — 1. The 
maximum of the specific heat (right panel) and its position 
T* (left panel) are reported as a function of h — h c where h c is 
the critical value of the external magnetic field h at the phase 
transition point. 

We report in Figj28]a detailed analysis of position and 
height of the low temperature peaks in the specific heat 
for all possible phase transition occurring at J = — 1. As 
already pointed out, the position T* of low temperature 



peaks is a linear function of the temperature while, con- 
trarily to what has been observed for the J = 1 case, the 
height remains constant. Furthermore, it is worth noting 
that, at the transition point, there is a jump in the height 
of the low temperature peak which can be traced back 
to the degeneracy ratio between the ground and the first 
excited states. As reported in Fig. this last finding, 
that represents a common fingerprint for all the phase 
transitions analyzed so far, is not observed at the tricrit- 
ical points because the first and the second excited states 
become degenerate with respect to the ground one in this 
very special case. 
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Figure 29: (Color online) Low temperature specific heat fea- 
tures at P2 tricritical point for n = 1 and J = — 1. Contour- 
plot and 3Dplot on the left report specific heat at as a function 
of h and T for U = 0. The four plots on the right report po- 
sition and intensity of low temperature peaks moving to the 
phase transition point fixing h and changing U (top) and vice 
versa (bottom). 

Entropy. The entropy can be calculated by means of 
two equivalent formulas: 



S(T, n) = S(0, n) + ^AdT 1 
S(T,n) = -f n ^ldn' . 



(51) 



The first expression requires the knowledge of the entropy 
at zero temperature which can be computed by means of 
the formula S(0,n) — fcgln(fi), where is the Boltz- 
mann's constant and f2 is the degeneracy of the ground 
state. Discarding few simple cases, it is generally not 
easy to compute fl and, in the case when the degeneracy 
of the ground state is not known a priori, the only way 
to calculate the entropy is the use of the second equation 
in (ISTj) . The study of the entropy plays a crucial role 
in the characterization and the identification of possible 
ordered phases in which we expect to have S — > in the 
thermodynamic limit. Filling and temperature depen- 
dence of the entropy for both J = 1 and J = — 1 cases 
are reported in Figs. I30ll32l 

For J = 1, as shown in Fig. [30] for different values of 
the filling, a finite ground state degeneracy is observed 
only in the NM-phase, as signalled by the value of the 
entropy in the T — » limit. In the F-phase, because of 
the presence of a spin ordering at T — 0, the entropy goes 
to zero at any filling as one would expect. Importantly, 
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as shown in Fig|3TJ when U is close to U c , S rapidly 
decreases showing a zero-temperature jump in the limit 
of f/ ->• U c . 




Figure 30: (Color online) Temperature dependence of the en- 
tropy plotted for different values of the external parameters 
corresponding to different phases at T — 0. 




Figure 31: (Color online) The entropy S at J = 1, h = 0.5 
as a function of the filling and for different values of T (left 
panel) and at n = 1 as a function of T for different values of 
U (right panel). 

For the case J = — 1, the filling dependence of the 
entropy is shown in Fig. [32] for different values of the 
external parameters U, h and T, corresponding to the 
phases observed at T — 0. As the temperature increases, 
the properties of the system at finite T can be described 
as a superposition of an increasing number of different 
configurations. In contrast, at low temperatures the en- 
tropy as a function of filling decreases rapidly near those 
values of n associated with charge and/or spin orderings. 
At these points, in fact, the state of the system can be 
uniquely described by a finite number of configurations 
(see FigsH ME and M for F, AF, Fl and F2-phase re- 
spectively) whose contribution to the entropy is expected 
to vanish in the thermodynamic limit. 

V. CONCLUDING REMARKS 

In this paper, we have evidenced how the use of the 
Composite Operator Method, based on Green's function 
and equations of motion formalism, leads to the exact 
solution of the t-U-J model in the narrow-band limit. In 
this limit the system can be described by a closed set of 
composite eigenoperators allowing for an exact solution 
that holds in principle for one, two and three dimensions. 
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Figure 32: (Color online) Entropy as a function of the filling 
for different temperatures and four sets of the external param- 
eters U and h corresponding to NM (top-left), AF (top-right), 
Fl (bottom-right) and F2 (bottom-left) phases observed at 
J = -1 and T = 0. 



We have considered the cases of both ferro (J > 0) and 
antiferro ( J < 0) couplings with an external longitudinal 
magnetic field. We have shown that a complete rigorous 
solution of the model does exist in the one-dimensional 
case in which Green's and correlation functions can be 
expressed in terms of four parameters, determined self- 
consistently. Within this framework, we have reported a 
detailed analysis of zero and finite temperature proper- 
ties for both ferro and antiferro exchange coupling J. We 
have found that in the T — > limit the system exhibits 
a very rich phase diagram characterized by long-range 
charge and spin orderings. In addition to the standard 
ferromagnetic and anti-ferromagnetic phases, that occur 
for J = 1 and J = — 1, respectively, in the case of antifer- 
romagnetic coupling the competition between J, U and h 
results in the presence of two ferromagnetic phases. We 
have reported a detailed analysis of each phase in terms 
of all the relevant single-particle correlators (magneti- 
zation, double occupancy, density-density and spin-spin 
CFs), thermodynamic quantities (free energy, entropy) 
and density of states. Furthermore, we have also re- 
ported that the presence of zero-temperature long-range 
orderings clearly characterizes the temperature depen- 
dence of response functions such as specific heat, charge 
and spin susceptibilities, leading to the presence of low- 
temperature charge and/or spin excitations. 

Summarizing, we have presented an exact and com- 
prehensive study of the one-dimensional t-U-J-h model, 
in the narrow-band limit. Our results show that this 
study is of substantial interest from many points of view. 
The model itself can be considered as a key-model for 
the study of phase transitions between different charge 
and spin orderings, generated by the competition among 
the different energy scales (U, J and h). The present 
exact solution is an important test for various approxi- 
mate approaches. Our exact solution of the model can 
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be used as a starting point for a perturbative expansion 
in powers of the hopping term in which the exact analyt- 
ical knowledge of the Green's function will be of crucial 
importance. This work is currently in progress. Further- 
more, we also remark the importance of the model from 
the point of view of statistical mechanics, motivated by 
the presence of several phase transitions, bicritical and 
tricritical points and anomalous behaviors in all the sys- 
tem response functions near the phase boundaries. More- 
over, because of the isomorphism discussed in Sec lIIDI 
for the one-dimensional case, the analytical framework 
developed in Sec |TT] also provides the exact solution for a 
two-leg ladder spin system with effective inter-chain and 
intra-chain spin-spin interactions in the presence of an 
external magnetic field. 
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Table I: Values of the coefficients A 
z = 2,4. 
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where a = (t,l) or a = (+1, —1) and I 5x5 is the 5x5 
identity matrix. The eigenvalues of these matrices are: 



Appendix A: Algebraic relations, spectral densities, 
energy and normalization matrices 

1. A^ coefficients 

By means of the algebraic relations ([5]) it is straight- 
forward to show that the filed n"(i) — - Yll=i n 3(*fc)> 
where if. (k = 1, . . . , z) are the first neighbor site of i, 
satisfies the following recursion formula: 
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2p-l 
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2d 



(Al) 
(A2) 



where are rational numbers, satisfying the sum rule: 
E A ^ = 1 • (A3) 



m— 1 

dp) _ 
,(p) 



For 1 ^ p ^ to, Am — S pm . For p > m the expressions of 



the coefficients Am depend on the coordination number 
z = 2d. For z = 2 we have: 
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In Tab.(|TJ we report some values for z — 2 and z = 4. 



2. Energy matrices 
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3. The spectral density matrices 



(s) 

The matrices fl a are defined as the matrices of rank 
(2z + 1) x (2z + 1) whose columns are the eigenvectors 

(s) 

of the energy matrices e a ■ It is worth noting that the 
fiS do not depend on a and s. For z = 2, from (|A6j) 
and (IA7I) we have: 
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The spectral density matrices pa'"^ are calculated by 
means of (TT51) in terms of the matrices ' and ijfi . The 
relevant elements of p^'^ for z = 2 have the expressions: 
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Pct:1,1 



(*,5) 
fV:l,l 



■51. 



4/ 



(^) 



^(7:1,3 + 4i <r;l,4 



4/ 



1 J ct;1,: 



4- 2T {s) 

T ZJ cr;l,3 



r(<0 



2/ 



00 



o-;l,2 



+ 2/ 



W 

ct;1,3 



j ct:1,4 zl t?;l,5 



As) 



i ( r;l,2 i o-:1.3 



•4/ 
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(All) 
(A12) 

(A13) 

,(A14) 

(A15) 
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4. The normalization matrix 

We recall that the normalization matrix is defined as 
Ia ] = ({^£°(i,*)> Vi s)t (i,*)})- By considering that this 
matrix is Hcrmitian and by recalling the recursion rule 
((6j) , it is easy to see that all the matrix elements can be 
expressed in terms of the elements of the first raw. By 
means of the algebra, these elements have the expres- 
sions: 

= K (p) _ a (p) 

Avf _ x(p) ^ (p = l,...,2z+l), (A16) 



<x;X,p 



where: 



«W = ([n%(i)r l ) , W = (n„ KMr 1 ) • (A17) 
For z = 2 the explicit expression of I a is: 



/ r(a) rW r(a) r« 



S J r(s) r(«) 
;2,1 J cr;2,2 J 



(X;l,3 i cr;i,4 



■(«) 



■W 

<r;l,5 
j(») 



rW t( s ) t( s ) t( s ) 



er;5,l • 



cr;4,l J cr;4,2 



<x;2,3 

<») , 

cr;3,3 - t <r;3,4 J cr;3,5 

r(«) 7-( s ) r( s ) 



<r;2,4 -V;25 
j(a) 



(a 



cr;4,3 

r« r(a) 

cr;5,2 J cr;5,3 J cr;5,4 



cr;5,4 ^er;4,5 

t;5,5 / 



(A18) 



with: 



<r;2,5 ~~ 
r(a) _ 
V;3,5 ~~ 
r(») _ 
V;4,5 ~~ 



(a) 



<r;5,5 



4 cr;l,2 ^ 
.lr(a) , 
4 J <t;1,3 ^ 
5 r(a!) 

" i6-V;i,2 

5 r(a) 
" 16 i o-;l,3 



5r(a) 

4 <t:1,4 ' 

5 r(a) 

4 cr;l,5 ' 
_ 21 r(a) 

16 <r;l,4 
_ 21 r(a) 
16 <7")1)5 



(A19) 



Appendix B: _Ho-representation 

By considering a one-dimensional system, let us write 
the Hamiltonian @ as: 



H = H, 



H 



flf - -2Jn 3 (i)n%(i) , 



I ' 



(Bl) 



where i is a site of the infinite chain. Since and H^p 
commute, for any operator O we can write: 



(0) = 



(B2) 



where the notation (. . . ). denotes the thermal average 



with respect tot he reduce Hamiltonian 



W. 



<•••>< = 



(B3) 



By using the algebraic properties (0 we have: 
e~ pH ^ = 1 + an 3 (*)fl(p) + [n(i) - 2D{i)} Q(p) , (B4) 



where: 

R(p) = n 3 (i 1 )+n 3 (i 2 ) + bn 3 (i 2 )[n{h)-2D(i 1 )] 



Q(P) 



with: 



+bn 3 (h) [n(i 2 ) - 2D{i 2 )} , 
6 [n(ii) - 2D{h) + n(i 2 ) - 2D{i 2 )] + 
+a 2 n 3 (H)n 3 (i 2 ) + 6 2 [n(h) - 2D(ti)] 
+ [n(i 2 ) - 2D{i 2 )] , 



a = sinh (J3J) 
b = cosh (13 J) 



1 



(B5) 



(B6) 



(B7) 



*p (p = 1, 2) are the first neighbor site of i. Now, we 
observe that H describes a system where the original 
lattice is divided in two non interacting sub-lattices (the 
chains to the left and right of the site i). Then, the 
correlation functions which relates two sites belonging to 
different sub-lattices can be written as the product of the 
two correlation functions: 



(a{j)b(m)) Q = (a(J)) (b(m)) 



(B8) 



for j and m belonging to different sub-lattices. By using 
these properties and the algebraic relations ([5]) we have: 

(e"^') = 1 + 2aG 2 A 2 [1 + b{X x - 2X 3 )} + 

+ (Gi - 2G 3 ) [a 2 A 2 + 2b {X x - 2X 3 ) + 
+b 2 (X 1 - 2X 3 f] , (B9) 

where we define: 



Gi = (n(i)>o 
G 2 = (n 3 (i)) 
G 3 = (D(i)) 



X, = (n(i + 1)) = <n(z - 1)) 
X 2 = (n 3 (i + l)) = (n 3 (i-l)) . 
X s = (D(i + l)) = (D(i-l)) 

(BIO) 

In order to calculate the quantities (n(i)) i , {n 3 {i)) v 
(D(i)) i let us consider the retarded Green's functions: 



G a £ '°\t-t') = (R[Ui,mUht')]) i 



G ( ^(t-t>) = (R[r 1a (i,tH(i,t>)]) i 



(Bll) 



In the iJo-representation the Hubbard operators satisfy 
the equations of motion: 



'imt® = KCO, «o(i)] = - (M + ha 3 ) 



iwAi) = Hi), H (i)] = - (M + ha 3 - U) r?(i) 



.9* 

By means of these equations we obtain: 

l-(ng(i)> 



gPV) 

Gi"^ 0) ( w ) 



aj+^-t-crh-f-i^ 
("^('))o 



u)-\-fj,-\-(7h — U +i8 

and for the equal-time correlation functions: 

l-(n 9 (i)) a 

G^' 0) - 



Co" 1 * ^ — v»0-(*)£ct(O)o _ I- 1 -, " : 

<^(0^(0) = i- (m ' (0>n 



(B12) 



(B13) 



(B14) 



19 



Recalling the algebraic relations: 

77^77* = Tt — (j - rifni , 
we obtain from (|B14[1 : 

- 13 " ( 1 + e 2 ,3 h + e 13 < M + ' l - 17 > ) 



G 2 



(«3(*)> i 



e /3h +e ,3n+ e /3(f + 2fc)+ e /3(2n + fi-<7) 
„f!(f I + 2fe)_ S M 



e fJh +e fln +e ^(n + 2h) +e f;(2 A i + fi-l7) 

r = tncw — e fs(2»+h-u) 



(B15) 



(B16) 



The parameters and can be straightforwardly 
calculated. For example: 



For a homogeneous solution, the internal parameters X\ , 
X 2 , X 3 are determined by means of the self-consistent 
equations: 



(n(i)> = (n(i p )) 
(n 3 (i)) = (n 3 (i P )) 
[(D(i)) = (D(i p )) 



(BIT) 



By using (|B4|) and the property (|B8|) . it is easy to show 
that: 



= Gi + 2aG 2 X 2 [1 + b(Xi- 2X 3 )] + 
+ (Gi - 2G 3 ) ■ [a 2 X 2 2 + 26 ■ 
■ (Xi - 2X 3 ) + b 2 (Xi - 2X3 ) 2 1 (B18) 



ra 3 (i)e 



i,(ii)e-^ w ) = Xi + aG 2 X 2 [l+Xi +26 (Xi- 2X3)] + 
+ (G 1 -2G 3 )[a 2 X 2 2 +6(l + X 1 )- 
■(Xi -2X 3 ) + 6 2 (X! -2X 3 ) 2 ] (B19) 

G 2 + 2a (Gi - 2G 3 ) X 2 [l + b (Xi 
-2X3)] + G 2 [a 2 X 2 2 + 26 (Xi - 2X 3 ) 
+6 2 (Xi - 2X 3 ) 2 ] (B20) 

(nsCiiJe-^^) = X 2 + (Gi - 2G 3 ) X 2 [ft + (a 2 + 6 + 6 2 ) 
■ (Xj - 2X3)] + aG 2 [xf (1 + 6) + 
+ (Xi - 2X3) + 6 (Xi - 2X 3 ) 2 ] (B21) 
G 3 (B22) 



D(i)e 



(1) _ 



K 

A« 



1 



^— {(Gi+aG 2 ) + 2a 



(B27) 



G 2 



2{e-P H Y 

+a (Gi - 2G 3 )J X 2 [1 + 6 (Xi - 2X 3 )] + 
+ (Gi - 2G 3 + uG 2 ) 



a 2 X 2 



-2b (X x - 2X3) + b 2 (Xi - 2X 3 ) 2 ] } , (B28) 



l^— {x 2 + aG 2 {(l + b)X 2 2 



+ [l + b(X 1 -2X 3 )} (X 1 -2X 3 )} + 
+ (Gi - 2G 3 ) X 2 \b + 

+ (a 2 + b + b 2 ) (X a - 2X3)] } , (B29) 



a! 2 ) = 



1 



2(e-f )H i 



(*) 



-{(Gi+aG 2 )X 2 



+a [G 2 + a (Gi - 2G 3 )] { (1 + b) X 2 + 
+ [l + b(X 1 -2X 3 )} (X 1 -2X 3 )} + 

+ (Gi - 2G 3 + crG 2 ) \bX 2 + 



(a 2 + b + b 2 ) X 2 (Xx - 2X3) 



(B30) 



k& = 



-L_ { Xx - 2X3 + X 2 + G 2 (/1 + 2/3) 



^(i + ^e-^W) = X 3 [l + aG 2 X 2 + 6 (Gi - 2G 3 ) 

(Xi - 2X3)] . (B23) 
Putting (|B18|) - (|B23|) into (|B17|) we obtain the equations: 

F x = G 1 -X 1 +aG 2 X 2 (l-X 1 ) + 6(Gi-2G 3 ) 

(1-X 1 )(X 1 -2X 3 )=0, (B24) 

F 2 — G 2 — X 2 + (Gi — 2Gs)X 2 

{2a - b + - [b + (a - b) 2 ] (Xi - 2X3)} + 

+G 2 {aX 2 (a - 6- 1) + (Xx - 2X 3 ) 

[26-a + 6(6-a)(Xi-2X 3 )]}=0, (B25) 

F 3 = G 3 -X 3 [l + aG 2 X 2 + b{G l -2G 3 ) 

{X 1 - 2X3)] = . (B26) 



+ (G x -2G 3 )(/ 2 + /4)} 



(B31) 



a! 3 ) = 



i^-^j (Gr + aG 2 ) {X 1 - 2X3 + X 



4 (e-^i 
+ [G 2 + a (Gr - 2G 3 )] (/1 + 2/3) + 
+ (Gi -2G 3 + aG 2 )(/ 2 + / 4 )} . 



(B32) 
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